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ABSTRACT 

A possible mechanism for screening of the surface magnetic field of an accreting neu- 
tron star, by the accreted material, is investigated. We model the material flow in the 
surface layers of the star by an assumed 2-D velocity field satisfying all the physical 
requirements. Using this model velocity we find that, in the absence of magnetic buoy- 
ancy, the surface field is screened (submergence of the field by advection) within the 
time scale of material flow of the top layers. On the other hand, if magnetic buoyancy 
is present, the screening happens over a time scale which is characteristic of the slower 
flow of the deeper (hence, denser) layers. For accreting neutron stars, this longer time 
scale turns out to be about 10 5 years, which is of the similar order as the accretion 
time scale of most massive X-ray binaries. 
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1 INTRODUCTION 

Radio pulsars can be broadly classified into two categories: 
: (a) isolated pulsars with rotation periods usually above 1 
s and very strong magnetic fields (10 12 - 10 13 G); and (b) 
binary /millisecond pulsars with much shorter rotation peri- 
ods and considerably weaker magnetic fields (10 s - 10 10 G). 
Observational indications have led to the proposal of a uni- 
fied model in which these two types of pulsars are regarded 
as two phases in the evolutionary history of neutron stars 
(see Bhattacharya 1995, 1996; van den Heuvel 1995; Verbunt 
& van den Heuvel 1995 and references therein). According 
to this model, a normal pulsar, if it is in a binary system, 
may undergo mass accretion from its binary companion for 
a limited period, and the angular momentum of the accreted 
matter may spin up the neutron star, ultimately making it 
appear as a millisecond pulsar after the accretion process 
stops. In this scenario, neutron stars seen in the X-ray as 
members of binary systems should correspond to the inter- 
mediate phase of transition from the normal pulsar phase to 
the millisecond pulsar phase. The fact that many millisec- 
ond pulsars are found in binary systems lends credence to 
this picture. If this scenario is correct, then the most natu- 
ral explanation of the weaker magnetic fields of millisecond 
pulsars must be that the magnetic field of a neutron star 
is significantly reduced during the accretion phase. Though, 
how this happens is still not understood very well. 



A variety of models have been proposed for the evolution of 
the magnetic field in an accreting neutron star (see Bhat- 
tacharya 1999a, 1999b; Konar & Bhattacharya 1999 and ref- 
erences therein) . There are two classes of models which have 
been mainly explored in this context: one that relates the 
magnetic field evolution to the spin evolution of the star 
and the other attributing the field evolution to direct effects 
of mass accretion. The different classes of models usually 
assume different kinds of initial field configurations. Mod- 
els depending on spin-down assume a core-flux supported 
by proton superconductor flux tubes. On the other hand, 
models invoking ohmic dissipation usually assume an initial 
crustal configuration. The mechanism of ohmic decay, be- 
ing unique to the crustal currents, is also used in models 
where spin-down is invoked for flux expulsion, for a subse- 
quent dissipation of such flux in the crust (Jahan Miri & 
Bhattacharya 1994; Bhattacharya & Datta 1996, Konenkov 
& Geppert 200, 2001a, 2001b). In an accretion- heated crust, 
the decay takes place principally as a result of rapid dissipa- 
tion of currents due to the decrease in the electrical conduc- 
tivity and hence a reduction in the ohmic dissipation time 
scale (Geppert & Urpin 1994; Urpin & Geppert 1995; Urpin 
& Geppert 1996; Konar & Bhattacharya 1997). 

However, in this context one of the least investigated mech- 
anisms that has often been invoked is the possible screen- 
ing of the field by accreted material. The idea of a possi- 
ble screening of the magnetic field of a neutron star by ac- 
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creting material was first proposed by Bisnovatyi-Kogan & 
Komberg (1974). Later, Taam & van den Heuvel (1980) indi- 
cated that the accreted matter, which is completely ionized 
plasma and hence diamagnetic in nature, might screen the 
pre-existing field, reducing its strength at the surface. Bland- 
ford, De Campli & Konigl (1979) suggested that the material 
accreting onto the poles will be confined by the strong mag- 
netic stresses near the surface of the star. At low accretion 
rates, the material sinks below the stellar surface until the 
hydrostatic pressure of the stellar material is as large as the 
magnetic pressure. The plasma then flows sideways giving 
rise to horizontal components at the expense of the verti- 
cal field and may result in a decrease in the observed field 
strength. Then it was shown by Woosley & Wallace (1982) 
and Hameury, Bonazzola, Heyvaerts & Lasota (1983) that 
the accretion column is like a small mountain on the polar 
cap rather than being subsurface. Semi-quantitative calcu- 
lations by Romani (1990, 1995) showed that hydrodynamic 
flows in the surface layers may bury the field to deeper lay- 
ers, effectively reducing the surface strength. 

When the accreting material flows horizontally from the 
polar caps to lower latitudes, magnetic field lines are ex- 
pected to be dragged by the flow, giving rise to distortions 
in the original field structure. This may lead to creation 
of additional horizontal components of the field at the ex- 
pense of the vertical ones. Now, horizontal magnetic fields 
are known to be subject to magnetic buoyancy (see, for ex- 
ample, Parker 1979, §8.7-8). Although the flow of accreting 
matter tries to screen the magnetic field, the field may pop 
back to the surface due to magnetic buoyancy (Spruit 1983) 
and Ohmic re-diffusion. The magnetic loops coming out to 
the surface due to buoyancy would reconnect with overlying 
field lines and thereby would make the screening ineffective. 
Simple estimates by Konar (1997) showed that the instabil- 
ity or the reconnection time scales happen to be too small 
for the screening mechanism to be effective. Similar con- 
clusion, as for the re-diffusion time-scale being of the order 
of 10 3 - 10 4 years, have been reached by Geppert, Page & 
Zannias (1999), in the context of the hyper-critical accretion 
onto a newly born neutron star and the re-emergence of its 
surface field. 

It seems at the first sight that the accreting material will 
not be able to screen the magnetic field of the neutron star. 
In particular, for higher rates of accretion, effects such as 
magnetic buoyancy and Ohmic diffusion might make the 
magnetic field emerge to the surface, rapidly, through the 
accreting matter. We, however, argue here that this is a 
complex problem of competing time scales and one has to 
do more detailed calculations, than what has been done pre- 
viously, to settle the question of screening of the magnetic 
field by accreted matter. To the best of our knowledge, no 
attempt has so far been made to investigate this problem in 
a 2-dimensional framework and ours is the first 2-D calcu- 
lation of the problem. Some of the results we present here 
depend crucially on the 2-D geometry and could not have 
been obtained by simpler 1-D calculations. 

The accreted material flowing from the polar caps would 
meet near the equator. This is expected to cause a submer- 
gence of such material under the surface over a belt around 



the equator. Even though it is difficult to say anything def- 
inite about the flow of the accreted material in the surface 
layers of the neutron star, without considering the Navier- 
Stokes equation including a reasonable viscosity; we do not 
really expect a pile-up of material in the equatorial region. 
Hence we expect that the submerging material would push 
against the solid interior of the neutron star and would dis- 
place the solid layers. Our rough estimates suggest the speed 
of this displacement of the solid interior to be ~ 10~ 6 cm s _1 
for Eddington accretion rates — comparable to the speeds 
of tectonic plates on the Earth's surface which are again 
~ 10~ 6 cm s _1 (see, for example, Emiliani 1992, §12.2). 
We propose that the solid inner material in the equatorial 
region, while remaining solid, is displaced very slowly down- 
wards as well as sideways to higher latitudes in the form 
of a counter-flow to the equator-ward flow of accreted ma- 
terial at the surface. This very slow movement of the solid 
material would carry the magnetic field with it, since mag- 
netic buoyancy would be ineffective in a region where the 
material is solid. The time taken by this very slow motions 
to cover an appreciable distance inside the neutron star is 
~ 10 5 years. Our 2-D calculations show that, in presence 
of accretion, the magnetic field can be screened in the time 
scale of this very slow interior motions of the solid material. 
Once the accretion is turned off, the submerged field should 
start re-diffusing to the outer layers. 

Although our present work is the first quantitative study of 
the effect on magnetic field of a flow in the meridional plane 
in the context of neutron stars, the effect of a meridional flow 
on solar magnetic fields has been studied in some detail (Dik- 
pati & Choudhuri 1994, 1995; Choudhuri & Dikpati 1999). 
Some recent models of the solar dynamo depend on the ex- 
istence of a meridional flow in the solar convection zone in 
a very crucial way (Choudhuri, Schiissler & Dikpati 1995; 
Durney 1995, 1997; Dikpati & Charbonneau 1999; Nandy 
& Choudhuri 2001). Our work derives its methodology from 
the models developed by the authors mentioned above. 

The mathematical formulation of the problem is discussed 
in the next section. Although we have run our code over a 
wide range of parameters, our code cannot handle the ex- 
treme values of parameters appropriate for the crust of a 
neutron star. The density varies in the outer layers of the 
neutron star by 8 orders of magnitude; the accreted mate- 
rial flows through a layer of thickness equal to 1% of the 
radius of the neutron star and the time scales of motion and 
Ohmic decay again differ by 8 orders of magnitude! In §3, 
we present some general results obtained with more moder- 
ate values of the parameters to understand the basic physics 
and show that one can draw some fairly generic, parameter- 
independent conclusions. The lessons learnt from this gen- 
eral study are extrapolated to the situation of the neutron 
star in §4. Finally, our conclusions are summarized in the 
last section. 

2 MATHEMATICAL FORMULATION 

The magnetic field of the neutron star will evolve according 
to the well-known induction equation: 
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^ = V x (v x B) - ^-V x (iv x B) 



(1) 



where a is the electrical conductivity of the medium (see, 
for example, Parker 1979; Choudhuri 1998). We assume an 
axisymmetric poloidal field which enables us to represent 
the magnetic field in the following form, 



B = V x [A(r,e)e^. 



(2) 



It is easy to show that the magnetic field lines are given 
by contours of constant r sin 9 A. On substituting eq.(0) in 
eq.(|lh, we find that A evolves according to the equation: 



OA 
dt 



+ I( v .V)(M) = r,(v 2 -^)A ; 



(3) 



where rj = c 2 /47t<t and s = rsin#. It is clear from eq.(^|) 
that only the poloidal component of v affects the evolution 
of A. 

To study the evolution of the magnetic field by solving 
eq.(^), we need to know the velocity field v. It is an ex- 
tremely non-trivial problem to calculate the velocity field v 
of the accreted material from the basic principles of fluid 
mechanics. We, therefore, follow the kinematic approach of 
specifying a velocity field v which has the required charac- 
teristics and then study the effect of this velocity field on 
the initial magnetic field. The material accumulated in the 
polar regions would cause an equator-ward flow at the sur- 
face in both the hemispheres. Let us assume that this flow 
is confined in a shell from r = r m to the surface r — r s 
(r m < r s), i- e i primarily in the liquid part of the crust. Near 
the equator, the flow originating from the two polar regions 
would meet, turn around and sink under the surface. We 
expect a pole- ward counter-flow in a region immediately be- 
low r = r m . Eventually this material settles radially on the 
core. So we expect v to become radially inward at some ra- 
dius r = rb (rt < r m ). The pole- ward counter-flow as well 
as the radially inward flow would take place mainly in the 
solid crystalline region. We have a source of new material in 
the polar cap region. Apart from that, we must have 



V.(pv) = 



(4) 



everywhere else. We now write down a velocity field which 
has all these characteristics. In the top layer r m < r < r s of 
equator-ward flow, we take 



pv r = K\ 
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In the next lower layer r b < r < r m where we expect the 
flow to turn around, we take 



pv r = Ki e 



-/3cos z 9 
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pv g = 




(8) 



Finally, when r < r b , we take the velocity to be radially in- 
ward (characteristic of the radial compression in the deeper 
layers): 



pv r 



pvg 



Ks 



0. 



(9) 



(10) 



It is evident from the context that the rate of accretion is 
related to these coefficients by the relation K3 = M /4-7T. The 
coefficients K\ , K2 and K3 can be related from the fact that 
pv r has to be continuous across r — r m and r = r b . From 
the continuity of pv r at r — r m , we obtain 



Ki 



= K 2 
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(11) 



whereas the continuity of pv r at r — rb gives 



I<2 

2 



-^r b + -r n 



+ 



K, 



= (I 2 ) 



We have specified pvg in such a fashion that it is equator- 
ward in the upper layer r m < r < r s and pole-ward in the 
lower layer r b < r < r m , falling to zero both at r — r m 
and r = r b . It is straightforward to show that V.(pv) = 
everywhere below r = r m . If we choose the values of (3 
and 7 in such a fashion that the factors er^yTJcos^) and 

^1 — e~ lB j have values different from unity in the non- 
overlapping regions around the equator and the pole respec- 
tively, then the divergence in the uppermost layer above 
r — r m turns out to be : 



V.(pv) 



p \ r m J sin 6* 



(13) 



We thus have a source of material only in the region around 
the pole in the upper layer. Our expression for the velocity 
field given by (p|)-(|lo"|) may seem somewhat complicated to 
the reader. But this is the simplest velocity field satisfying 
all the necessary characteristics that we have been able to 
write down. 

Fig. jl] shows how pv, given by eqs.(P)-(p"o|), looks like. In all 
our discussions, we shall assume the radius r s of the neu- 
tron star to be the unit of length. Fig. |l| has been generated 
by taking r m = 0.75, r b = 0.5, f3 = 1.0 and 7 = 1.0. If 
we assume the density to be constant, then the amplitude 
of the equator-ward flow in the top layer and the ampli- 
tude of the pole-ward counter-flow in the layer underneath 
it are of comparable magnitude. But in the crustal layers of 
a neutron star, the density varies by 8 orders of magnitude, 
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Figure 1. The top panel shows pv in a right-angular slice (0 < 
6 < tt/2, 0.25 < r < 1.0). We have used r m = 0.75, r b = 0.5 in 
this picture. The bottom panel shows the contours of V • (pv) for 
this flow velocity. Note that the divergence vanishes everywhere 
except in a narrow region near the polar cap. 



which is difficult for a numerical code to handle. Therefore, 
to study the effect of density stratification, we allow the 
density to vary by 2 orders. This is done by modeling the 
density variation with the radius as 



p{r) 



1.01 



-a(r s — r) 



(14) 



with a = 5.0. A plot of this variation of density with radius 
is shown in Fig. This figure also shows the variation of vg 
as a function of r at the mid- latitude ($ = 7r/4). The unit of 
velocity is set by equating Kz appearing in eq.(^) to unity. 
This, along with the unit of length defined by the radius of 
the neutron star, fixes the unit of time in our problem. So, 
the unit of time is essentially the convective time scale from 
eq-(0), given by 



(15) 
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Figure 2. Radial dependence of density and horizontal velocity 
vg at mid-latitude. As we move deeper below the surface, note 
that vg becomes smaller, changes sign and then drops to zero 
as the flow velocity becomes increasingly radial in the interior 
regions. 



We see in Fig. || that the equator-ward flow velocity near the 
surface is about 100.0. The time scale for the flow (i.e. the 
time taken by this flow to traverse the radius of the neutron 
star) is about 0.01 unit of time. The pole-ward counter-flow 
under the top layer has an amplitude of about —2.0, giving 
a time scale of about 0.5. This flow is weaker by a factor 
of about 100 compared to the flow near the surface because 
the density in this deeper layer is enhanced by a factor of 
100 with respect to the density at the surface. 

To obtain the evolution of the magnetic field with time, 
we need to solve eq.(^) with the velocity field defined by 
eq.(p|)-(|ic|). We integrate eq.(^|) numerically in the region 
< 8 < 7r/2 and rrj < r < r s , where the interior boundary 
of integration ro is taken some distance below rt- The follow- 
ing boundary conditions are satisfied. The field lines from 
the two hemispheres should match smoothly at the equator, 
which requires 

r)A 

00=0, at 9 = tt/2. (16) 

To avoid a singularity at the pole, we need to have 

A = 0, at 61 = 0. (17) 

Now, at the upper boundary r = r s , we allow the magnetic 
field to smoothly match a potential field outside (Dikpati 
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& Choudhuri 1994), which satisfies the free space Laplace 
equation 

2 1 



0. 



Its solution has the general form 
A(r > r„,6,i) 



(18) 



(19) 



which, to match the boundary condition at 8 — n/2 men- 
tioned above, picks up only the odd n terms. The coefficients 
a n are obtained from the surface values of A by the following 
relation 



2n + l 
a„{t) = — — — —r B 



r/2 



n(n + l) JQ 

The derivative of A at the surface is 
dA 



A(r s , 6, i)P„ (cos 0) sin 6 d6 .(20) 



dr 



(21) 



Hence, the top boundary condition is that our solution has to 
match (21) at r = r s , with a n (t) being given by (20). Finally, 
at the lower boundary r — ro, we allow the magnetic field to 
be carried freely downward by the radially inward velocity 
field there. How this is implemented in the numerical scheme 
is pointed out in the Appendix. 

Our code is developed by using the code of Dikpati & 
Choudhuri (1994) developed in the context of solar MHD 
as a template. We summarize the salient features of the nu- 
merical schemes in the Appendix of the present paper. 



3 FIELD EVOLUTION - GENERAL RESULTS 

So far, the evolution of a magnetic field due to a flow like 
what we are considering has not been studied before. There- 
fore, to begin with, we present a general study of the prob- 
lem, to gain an understanding of it, before we consider ap- 
plications to the neutron star. In a neutron star, the veloc- 
ity field remains confined to the uppermost layer, which is 
rather thin compared to the size of the star (the equator- 
ward flow takes place in a layer of 100 m, which is about 
1% of the neutron star radius). Solving the equations nu- 
merically in such a thin layer involves dealing with special 
numerical instabilities. Nor is it easy to plot the resulting 
magnetic field configurations without making some trans- 
formations to blow up the thin region where all action takes 
place. We, therefore, first study the basic physics of the 
problem by taking the velocity field confined to a reason- 
ably thick layer. Most of the calculations presented in this 
section uses the velocity field shown in Fig. |l| We shall see 
that it is possible to draw important conclusions which are 
independent of the thickness of the flow region. The impli- 
cations of our calculations for an accreting neutron star will 
be discussed in the next section. 




Figure 3. Initial field configuration, assuming a central dipole 
(given by eq. (p^)), in a region defined by 0.25 < r < 1 and 
< 9 < 7I-/2. 




Figure 4. Field configuration at t = 0.1 starting from the initial 
configuration shown in Fig. 3 evolved with r) = 1. 



to be in the form of a dipolar potential field centred at the 
centre of the star and pervading the entire star. Fig. ^ shows 
the field lines for such a potential field, given by 



A = 



sin 6 



(22) 



which we assume to be the initial configuration of the field. 



3.1 A Central Dipole 

There is no clear consensus regarding the question as to 
whether the magnetic field of a typical neutron star pene- 
trates through the core or whether it is confined to the crust. 
We shall first present results taking the initial magnetic field 



3.1.1 Without Magnetic Buoyancy 

We have already mentioned that the magnetic fields in the 
uppermost molten layers just below the surface would be 
subject to magnetic buoyancy. Let us first present some cal- 
culations without including magnetic buoyancy. Now, to cal- 
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Figure 5. Field configuration at t = 0.1 starting from the initial 
configuration shown in Fig. 3 evolved with rj = 0.01. 
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Figure 6. Evolution of the mid-latitude surface field with time. 
The curves 1 and 2 correspond to the cases of Fig. ^ (r) = 1) and 
Fig. 5 (r) = 0.01) respectively. 
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culate the evolution of the field in time, we also need to 
specify the value of rj. Figs. ^ and |B| show the magnetic field 
configuration after time 0.1 for rj — 1 and rj = 10~ 2 respec- 
tively (in units in which r s and Kz are taken as unity and 
we have also suppressed the physical dimensions of rj here). 
When rj has the smaller value 10 -2 , the magnetic field is 
nearly frozen in the plasma and gets buried underneath the 
layer of accreted material, as seen in Fig. ||. However, if rj has 
the larger value of 1, then the field emerges out to the sur- 
face due to Ohmic diffusion and it is no longer possible for 
the accreting material to bury the magnetic field very effec- 
tively (Fig. W) . The time evolution of B r at a pointjust below 
the surface at the mid-latitude is shown in Fig. pi It is seen 
that B r decreases more rapidly if the diffusivity is low. At 
the first sight, this may appear counter-intuitive. We would 
normally expect a lower diffusivity to make a magnetic field 
decay less rapidly. Here the situation is reversed because the 
flow of accreting material is able to bury the magnetic field 
better if rj is lower and the magnetic field is frozen in the 
material more effectively. From a plot like Fig. [], one can 
easily estimate the time 7i/ 2 in which the magnetic field 
near the surface falls to half of its initial value. Figure 7 
shows how T1/2 varies with rj. When rj is large, the flow can- 
not carry the magnetic field with it and Ti/2 is large. On 
the other hand, a smaller rj makes the magnetic field frozen 
in the plasma and the field can be carried away in the time 
scale of the flow. Hence, for smaller values of rj, Tx/ 2 reaches 
an asymptotic value which is essentially the time scale of 
flow in the upper thin layer. The fact that Tx/ 2 reaches an 
asymptotic value even before the value of rj is reduced to 
10" 2 is of some significance. Our code becomes numerically 
unstable if we make rj an order smaller than 10 -2 . On the 
basis of Fig. [?], we expect that the behaviour of the system 
should not change on lowering the value of rj below 10 -2 , 
even though we are unable to run our code for such values. 

3.1.2 With Magnetic Buoyancy 

In the uppermost layers of a neutron star where matter is in 
a molten state, magnetic buoyancy would make the magnetic 
field rise against gravity. It has been shown in the context of 
the solar convection zone that, when rotation is important, 
the magnetic field rises parallel to the rotation axis rather 
than radially outward (Choudhuri & Oilman 1987, Choud- 
huri 1989). In a rapidly spinning neutron star, the magnetic 
field may tend to rise parallel to the rotation axis. But, to 
simplify the problem, we are assuming that the magnetic 
field rises radially. The results of the present calculation are 
not expected to change much if we make the rise parallel to 
the rotation axis. 

We capture the effect of magnetic buoyancy by superimpos- 
ing on our velocity field given by eq.(g-^o|) a radial velocity 
in the uppermost layer. Hence, in the top layer r m < r < r a , 
we add a velocity given by 

v r = v mb (r - r m ). (23) 



Figure 7. Variation of T 1/2 for the surface field with n when the The parameter v mh determines the strength of magnetic 

initial field configuration is given by Fig. B. buoyancy. It should be mentioned here that this is just a 

simplified way of modeling the effect of magnetic buoyancy. 
In reality, the upward velocity due to buoyancy must be 
proportional to the square of the field strength at any point 
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Figure 10. Evolution of the mid-latitude surface field with time. 
The curves 1 and 2 are for cases without and with magnetic buoy- 
ancy corresponding to Fig. |] and Fig. |8 respectively. For both the 
cases, r\ = 0.01 and u m b = 50 for curve 2. 



Figure 8. Field configuration at t = 0.1 starting from the ini- 
tial configuration shown in Fig. 3 evolved with rj = 0.01 in the 
presence of magnetic buoyancy. Here, we have taken u m b = 50. 




Figure 9. The field configuration of Fig. H evolved to t = 0.5. 



(Cumming, Zweibel & Bildsten 2001). Since, we are only 
interested in obtaining the gross features of the possible 
screening effect, here we work with the above mentioned 
model velocity instead of going into further details. We have 
found that our results remain virtually unaltered on taking 
u m b equal to 50.0 or 100.0. This means that the effect of 
magnetic buoyancy saturates when u m b = 50.0 and things 
do not change any more on increasing v m b. We now present 
calculations with « m b = 50.0. We again begin with the ini- 
tial configuration given by eq.Q and shown in Fig. g. Now, 
Fig. |shows the field at t = 0.1 for r? = 10" 2 . The only differ- 
ence between Fig. g and Fig. |^ is that Fig. |^ is obtained with 
magnetic buoyancy included. On comparing the two figures, 
we find that magnetic buoyancy in the uppermost layer 



has the effect of making magnetic field lines radial there. 
In Fig. ^, we have seen clearly that the equator-ward flow 
completely screens the magnetic field from the uppermost 
surface layer. This has not happened in Fig. |§| where mag- 
netic buoyancy has made the magnetic field emerge through 
the uppermost layer. To show how the magnetic field will 
keep evolving from the configuration of Fig. ^, we show the 
magnetic configuration at a later time 0.5 in Fig. ^. It is 
clear from Figures ^ and ^ that the slow flow in the deeper 
layers carries the magnetic field deeper down as well as to 
higher latitudes from the equatorial region. How the mag- 
netic field B r at mid-latitude near the surface changes with 
time is shown in Fig. [w| For the sake of comparison, in this 
figure we again show the evolution of B r without magnetic 
buoyancy for rj = 10~ 2 (which was shown in Fig. ^ also). 
Although the magnetic field in the case with buoyancy has 
not been screened by the equator-ward flow in the upper 
layer, as in the case without buoyancy, we still find that B r 
keeps decreasing even in the case with buoyancy. The time 
in which B r decreases appreciably is of the order 0.1. We 
have already estimated in §2 that the time scale for the slow 
flow in the interior is 0.5. This clearly indicates that the de- 
crease in B r near the surface is caused by the magnetic field 
being carried to deeper layers by the slow flow in the inte- 
rior. Even though magnetic buoyancy is now allowing the 
magnetic field to emerge through the top layer and it is not 
possible for the equator-ward flow in that layer to screen the 
magnetic field, we conclude that the magnetic field still gets 
screened from the outside by the effect of the slow interior 
flow which pushes the magnetic field deeper down. In other 
words, the magnetic field gets buried in the time scale of 
the slow interior flow when magnetic buoyancy is taken into 
account (instead of getting buried in the time scale of the 
faster flow of the top layer as in the case when buoyancy is 
not included). 

Letting the field evolve in time for different values of 77, we 
can find Ti/ 2 for the different values when magnetic buoy- 
ancy is included. Fig. [y] shows how Tx/ 2 varies with 77 in 
this case. As in Fig. M, we find that 7\/2 becomes indepen- 
dent of 77 when 77 is sufficiently small to make the magnetic 
field nearly frozen. The asymptotic value of T1/2 in Fig. |h| 
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Figure 11. Variation of Tj/ 2 for the surface field with 77 when 
the initial field configuration is given by Fig. |jj and in presence of 
magnetic buoyancy (i> m b = 50). 




Figure 13. Field configuration at time t = 0.1 starting from the 
initial configuration shown in Fig. ^ evolved with rj = 0.01 and 
without magnetic buoyancy. 




Figure 12. The initial configuration of the crustal field, in a 
region defined by 0.25 < r < 1 and < 9 < tt/2. 

however, is much larger than what it is in Fig. [jj. The asymp- 
totic value of Tyji in Fig. ^ basically corresponds to the time 
scale of the equator-ward flow in the top layer, whereas the 
asymptotic value in Fig. |ll] corresponds to the time scale of 
the slower flow in the interior. 




Figure 14. Field configuration at time t = 0.1 starting from the 
initial configuration shown in Fig. ^ evolved with r\ = 0.01 and 
with magnetic buoyancy (u m b = 50). 



3.2 Crustal Field Configuration 

We now show what happens if the magnetic field is confined 
to the crust. For this purpose, we begin with the magnetic 
configuration shown in Fig. Then Figures jl3| and 14| re- 
spectively show the field configurations at t = 0.1 without 
and with magnetic buoyancy (using 77 = 10 -2 ). It is again 
found that the flow of accreted material in the top layer 
quickly screens the magnetic field in the absence of mag- 
netic buoyancy, whereas on inclusion of magnetic buoyancy 
we find that the field is able to emerge through the top layer. 
The interior flow carries the magnetic field slowly to deeper 
layers and to higher latitudes. It is interesting to note that 



the magnetic field eventually tends to take up configura- 
tions which do not look very different from the configura- 
tions which arise if we begin with an initial magnetic field 
passing through the core (compare Figures Il3l— IL4l with Fig- 
ures |E| and |§]) . We thus expect the results to be qualitatively 
similar even if we start from quite different initial configu- 
rations. 

To see the effect of the velocity field very clearly, we also 
study what happens if the velocity field is switched off and 
the mag netic field is allowed to decay due to diffusion alone. 
Fig. n5l shows what the magnetic field starting from the ini- 
tial configuration of Fig. |l2j would look at t = 1.0 if the 
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netic field would have been screened in a very short time 
scale arising from the flow of accreted material in the top 
layer. On account of magnetic buoyancy, we obtain an inter- 
mediate time scale between these two, which is essentially 
the time scale of the slow interior flow. 




Figure 15. Field configuration at time t = 1.0 starting from the 
initial configuration shown in Fig. |l^, in the case of pure diffusion, 
with jj = 0.01. 




Figure 16. Evolution of the mid-latitude surface field with ti 
The curves 1, 2 and 3 correspond to the cases seen in Fig. |1 
Fig. [14] and Fig. |l5| respectively. 



velocity field is switched off. Since diffusion is a slower pro- 
cess for the values of parameters we are using, note that we 
have evolved the magnetic field for a longer interval com- 
pared to what was done in the other figures. We see that 
the magnetic field has spread outside the region where it 
was initially confined. We point out that such a study of 
pure diffusion would not be possible with the initial poten- 
tial field passing through the core (shown in Fig. ^ that we 
considered in the earlier subsections. Our bottom boundary 
condition specified in the Appendix is such that A would 
not change at the bottom if the velocity field were switched 
off. Hence, if we start from the initial potential field (given 
by eq.(^)) and allow it to evolve under diffusion alone, it is 
easy to see that the field remains invariant. 

Fig. |l^ shows how B r at mid-latitude near the surface 
evolves for the three cases shown in Figures |l3| - |l5| (i.e. with- 
out magnetic buoyancy, with magnetic buoyancy and with- 
out any flow). We see that pure diffusion has a larger time 
scale. If there were no magnetic buoyancy, then the mag- 



3.3 Field Confined to a Narrow Layer 

We now argue that our general conclusions do not change 
on reducing the thickness of the layer in which the flow 
takes place. We point out in the next section that we expect 
r-m, = 0.99 in a neutron star (we have been taking r m = 0.75 
in the above exploratory calculations) . We find that our code 
does not retain numerical stability if r m is increased to 0.99. 
However, we have been able to obtain results by increasing 
the value of r m up to 0.95 and find that our general conclu- 
sions inferred in the earlier subsections do not change. We 
now present our results obtained by taking r m ~ 0.95 and 
rt — 0.85. The initial field configuration is shown in Fig. [l7j 
We also change the value of a in eq.(|l^) to 20 to simulate a 
very sharp drop in the density profile. 

We find that in absence and in presence of magnetic buoy- 
ancy the field configuration gets modified to the ones seen in 
Fig. [l^ and Fig. ^ respectively. It is evident that the generic 
nature of field evolution is very much the same as that seen 
in the previous subsection. In the earlier case, however, mag- 
netic buoyancy saturated at D m b = 50 and the nature of field 
evolution did not change on increasing n m b beyond. In the 
present case, the similar effect is seen around w m b = 500. 
This difference is not difficult to understand if we keep in 
mind that « m b has to be multiplied by (r — r m ) to get the 
buoyant rise speed, as seen in (23). Whereas the maximum 
value of (r — r m ) was 0.25 for the calculations in the previous 
subsections, here its maximum value is 0.05. To compensate 
for this, we need to take a larger w m b to get buoyant rise 
speeds of similar magnitude. 

For comparison, we also show the case of pure diffusion in 
Fig. ^(j. To compare the time scales of evolution, in Fig. ^l] 
we have plotted the decay of the surface field (near mid- 
latitude) for the three cases - a) without magnetic buoyancy, 
b) with magnetic buoyancy and c) for pure diffusion. Con- 
forming to our expectation (which has already been seen 
in the previous sub-section), the time scales of field decay 
without buoyancy, with buoyancy and that of pure diffusion 
arrange themselves in an increasing order. Therefore, in con- 
clusion we can say that all the characteristic features of our 
calculation are also borne out by confining the physics to a 
very narrow layer near the surface. This vindicates the at- 
tempt of generalizing our results to the case of neutron stars 
(which we describe in some detail in the next section), even 
though the exact calculations could not be performed by our 
numerical code. 



4 FIELD EVOLUTION IN NEUTRON STARS 

In a real neutron star, the equator-ward flow must be con- 
fined to a very narrow layer near the surface within which 
density varies by several orders of magnitude. Our numeri- 
cal code is unable to handle such a severe density contrast. 
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Figure 17. Initial field configuration similar to that in Fig. 
but confined to a narrower region with r < 0.8. 



Figure 19. Field configuration at time t = 0.1 starting from the 
initial configuration shown in Fig. [l^ evolved with r) = 0.1 and 
with magnetic buoyancy (« m b = 500). 



Figure 18. Field configuration at time t = 0.1 starting from the 
initial configuration shown in Fig. [l^ evolved with rj = 0.1 and 
without magnetic buoyancy. 



Figure 20. Field configuration at time t = 0.5 starting from the 
initial configuration shown in Fig. |l^, in the case of pure diffusion, 
with rj = 0.1. 



However, from the results presented in the previous section, 
the following important parameter-independent conclusions 
emerge : 

(i) If the flow time scale is shorter than the diffusion time 
scale, then it is possible for the flow to carry the magnetic 
field with it and eventually bury the field under the surface. 

(ii) In the absence of magnetic buoyancy, the magnetic 
field gets buried in the time scale of the equator-ward flow 
in the top layer. 

(iii) When magnetic buoyancy is taken into account, the 
magnetic field is buried in the time scale of the much slower 
counter-flow under the top layer. 



These results seem fairly robust and we expect them to be 
valid in the case of neutron stars too. Therefore, we are able 
to make some inferences about the neutron stars on the ba- 
sis of our estimates of the time scale of diffusion, the time 
scale of equator-ward flow of accreted matter and the time 
scale of very slow motion of the solid material under the top 
layer. We now attempt to estimate these various time scales. 

The magnetic field of the isolated neutron stars is typically 
~ 10 12 G, which makes the magnetic pressure at the polar 
cap ~ 10 23 dyne cm -2 . As material starts collecting in the 
polar regions, excess pressure builds up, and when this ex- 
cess pressure becomes ~ 10 2 of the magnetic pressure, the 
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Figure 21. Evolution of the mid-latitude surface field with time. 
The curves 1, 2 and 3 correspond to the cases seen in Fig. [IsJ 
Fig. [l9| and Fig. [20] respectively. 



accumulated material flows sidewise overcoming the mag- 
netic stresses (Brown & Bildsten 1998). From the equation 
of state of the neutron star material, we find that the den- 
sity at which the pressure becomes ~ 10 25 dyne cm -2 is 
~ 10 9 gem -3 . From the density profile of a 1.4 M0neutron 
star, assuming a moderately stiff equation of state, we find 
that p — 10 9 gem -3 corresponds to a depth of 100 m below 
the surface, where r has the approximate value ~ 0.99r s (r 3 
is radius of the star, which turns out to be about 10 km). 
(For a discussion on the equation of state used here and the 
calculation of the density profile of a 1.4 Mq neutron star, 
see Konar 1997). We expect that the excess pressure of ac- 
creting material at the poles will induce equator-ward flow 
up to this depth of 100 m. 

Let us assume that the rate of material accretion in the polar 
region is M. If the polar cap has a circumference of 2nl and 
there is an outflow across this circumference with velocity v 
and up to depth h, then we must have: 



M fa 2-Klhpv . 



(24) 



The Eddington accretion rate for a neutron star is ~ 10 
M Q yr -1 , i.e. about 10 18 gm s _1 in cgs units. Taking I fa 1 
km, h fa 100 m, p fa 10 9 gm cm -3 , we find the velocity of 
outflow to be: 

v fa 10 _1 cm s~ x . (25) 

Combining this velocity with a length of 10 km (i.e. the typ- 
ical radius of a neutron star), we get a time scale of about 
1 year for the equator-ward flow of accreted matter. 

We now estimate the velocity of the slow counter-flow of 
solid material under the surface. Assuming a density increase 
of 100, we found that the velocity of counter-flow in the inner 
layers is about 100 times less than the velocity of equator- 
ward flow in the top layer (see Fig. || and discussions in the 
text about it). We have estimated the density in the layer 
of equator- ward flow to be about 10 9 gm cm -3 , whereas the 
density some distance below the surface is ~ 10 14 gcm~ 3 . 
An increase of density by a factor of 10 5 would lead to a de- 
crease of velocity by the same factor. Although geometrical 
factors such as the radial depth of the region across which 
the counter-flow takes place would have some role, here we 



are interested only in rough orders of magnitude. Dividing 
the velocity at the surface by this factor of 10 5 , we get a 
velocity of about 10 -6 cm s -1 . This corresponds to a time 
scale of ~ 10 5 years. We can conclude that the magnetic field 
of the accreting neutron star would get screened in 1 year 
in the absence of magnetic buoyancy and in about 10 5 years 
if magnetic buoyancy is important, provided the diffusion 
time scale is larger than these time scales and the magnetic 
field can be assumed as nearly frozen in the material. 

The electrical conductivity a in the interior of an accretion 
heated neutron star is ~ 10 26 s _1 so that r\ ~ 10~ 6 cm 2 s -1 . 
Using a length scale of about 1 km, we obtain a diffusion 
time scale ~ 10 16 s or ~ 10 9 years. It is true that the elec- 
trical conductivity of the crust can be as low as ~ 10 20 s _1 
in the surface layers and the diffusive time scale calculated 
on the basis of that will be smaller by several orders. How- 
ever, if magnetic buoyancy is dominant in the outer layers, 
the magnetic field comes out to the surface in any case and 
we need not bother about the role of Ohmic diffusion in 
the crust in making the magnetic field emerge through the 
accreted matter. The important point to note here is that 
the diffusion time scale in the interior is much larger than 
the estimated time scale of 10 5 years of the slow flow in the 
interior. We, therefore, expect that the slow displacement 
of the solid material in the interior to be able to carry the 
magnetic field and bury it in the time scale of 10 years. It 
should be mentioned here that this time scale happens to be 
of the similar order as the mass accretion time scale of most 
massive X-ray binaries. 

Let us now make a few comments on magnetic buoyancy. It 
is well known that magnetic buoyancy is particularly desta- 
bilizing in a region of super-adiabatic gradient, whereas its 
effect can be reduced in a region of sub-adiabatic gradient 
(Parker 1979, Moreno-Insertis 1983). Although the top layer 
of the neutron star is expected to be sub-adiabatic, we be- 
lieve that magnetic buoyancy will still be important there. 
In a sub-adiabatic region, a magnetic flux tube continues 
to rise as it exchanges heat with the surrounding medium 
and comes to thermal equilibrium (Parker 1979, §8.8). The 
very high thermal conductivity in the crust of the neutron 
star should keep magnetic buoyancy going even in a sub- 
adiabatic environment. So it is unlikely that the magnetic 
field of the neutron star can get screened in the time scale of 
1 year of the flow of accreting matter on the surface. Even 
if magnetic buoyancy were suppressed for some unknown 
reason, the much shorter Ohmic diffusion time in the crust 
would make the magnetic field emerge out of the accreted 
matter (Cumming, Zweibel & Bildsten 2001). However, it 
appears that the slow motion of the solid material in the in- 
terior of the neutron star may carry the magnetic field with 
it and eventually bury the magnetic field in its time scale of 
10 5 years. It is true that we have not been able to present 
detailed calculations appropriate for the neutron star and 
some of our arguments are heuristic. However, we humbly 
submit that this proposal is worth considering. 
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5 CONCLUSION 

In the past there has existed a controversy regarding the 
possible mechanism of screening of the magnetic field of the 
neutron star by the accreting matter. We have tried to ad- 
dress the question with a 2-D simulation. Although we could 
run our code only for thicknesses of flow layers much more 
than what they would be in a realistic situation, we have 
arrived at a few generic results which seem robust and in- 
dependent of the thickness of flow layer. So we believe that 
these results can be extrapolated to neutron stars. Basically 
we find that the magnetic field would get screened in the 
short time scale of the flow of accreting material in the top 
layer if magnetic buoyancy is neglected, whereas on inclu- 
sion of magnetic buoyancy the screening takes place in the 
somewhat longer time scale of the slow interior flows. Since 
magnetic buoyancy is expected to be important for neutron 
stars, we believe that the second time scale of interior flows 
will be the appropriate time scale for screening the magnetic 
field of neutron stars. We have estimated this time scale to 
be of order 10 5 years. 

It is remarkable that this is comparable to the time scale of 
accretion in massive X-ray binaries. It must be an extraor- 
dinary coincidence that these two time scales arising out 
of completely different considerations turn out to be fairly 
close. This coincidence has important implications. If the 
time scale for screening was much larger than the time scale 
for accretion, then the magnetic field could not be screened 
much during the accretion phase. On the other hand, if the 
screening time scale happened to be much shorter than ac- 
cretion time scale, then the magnetic field might have been 
screened so well that there would be virtually no magnetic 
field seen from outside. Observationally we find that the 
magnetic fields of binary /millisecond pulsars are consider- 
ably weaker compared to magnetic fields of normal pulsars, 
but are not altogether negligible. Can we take this as an ob- 
servational confirmation that the screening time scale should 
be of the similar order as the accretion time scale? 

Although considerations of neutron stars motivated us to 
undertake these calculations, we point out that what we 
have studied is the basic physics of magnetic field evolution 
in a situation of polar accretion. Our results should apply 
to other situations such as accretion onto white dwarfs. In 
the case of neutron stars, we have observations which seem 
to indicate that accreting matter may screen the magnetic 
field. We do not have observations to tell us what happens 
to magnetic fields in other stellar systems as a result of mass 
accretion. But the basic physics must be the same. 

Our model is based on certain assumptions. We have as- 
sumed a specific form of the velocity field, which essentially 
implies that the accreting matter primarily sinks near the 
equator where flows from the two different hemispheres meet 
and the matter in the underlying layers at the equator is 
thereby pushed downward as well as to higher latitudes in 
the form of a counter-flow. The other crucial assumption is 
that the dense matter in the lower layers is displaced very 
slowly while remaining solid and hence the magnetic field 
remains embedded in the solid matter of interior flow re- 
gions where magnetic buoyancy is suppressed. Once we make 



these assumptions, our results follow from fairly straightfor- 
ward calculations. Are our assumptions justified? We are 
not aware of any arguments that can be advanced against 
these assumptions and they do appear reasonable to us. To 
substantiate these assumptions fully, one will require some 
sophisticated fluid mechanics and sophisticated condensed 
matter physics of very dense material. We hope that re- 
searchers more conversant in these subjects than us will look 
into these questions in future. 

This being the first calculation of this type, we have re- 
stricted ourselves to a simple exploratory model. We could 
not run our code reducing the thickness of the upper layer 
to realistic values. Researchers more competent in handling 
numerical instabilities should be able to do better than us. 
Also, we have made the simplifying assumption that our ve- 
locity field remains independent of time. As magnetic field is 
pushed away, the polar cap widens and the accretion should 
take place over a broader region around the pole. We do not 
expect this to change our basic results, which seem quite ro- 
bust. However, one will have to use a time- varying velocity 
field to build more realistic models of diamagnetic screening. 
All our calculations are based on a velocity field for which we 
could write down a convenient analytical expression. Since 
we can only make intelligent guesses as to what the velocity 
field may be like at the surface of the neutron star, it would 
be worthwhile to investigate other forms of the velocity field 
and to ascertain whether certain velocity fields are more effi- 
cient in diamagnetic screening than others. We are right now 
in in the process of studying the screening of the magnetic 
field by other types of velocity fields, especially by veloc- 
ity fields which change with time as the polar cap widens. 
We hope to present our results in a future paper (Konar 
& Choudhuri, in preparation). We end with the comment 
that our first exploratory study appears promising and the 
scenario presented by us should be investigated further. 
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APPENDIX 

The basic features of our numerical code developed to solve 
eq.Q are outlined here. The basic equation we have to solve 
is (3), which can be written down in the form 

— = L r A + LeA, (26) 
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where the operators L r and Lg are given by 
L r A = 
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(28) 



To solve this two-dimensional problem numerically, we em- 
ploy the Alternating Direction Implicit (ADI) scheme of dif- 
ferencing (see, for example, Ames 1977; Press et al. 1988). 
In this scheme, each time-step is divided in two halves. In the 
first half time-step one direction (say, r) is advanced implic- 
itly and the other direction (say, 6) is advanced explicitly. In 
the next half time-step the two directions are treated in the 
reverse manner. If Afj is the value of A at the grid point 
(i, j) at time step m, then this implies the following two time 
steps: 
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(30) 



Here L r and Lg are the difference forms of the operators L r 
and Lg. They have the forms 

L r A™ 3 = a{i,j)AT-i,j +b(i,j)A™ j + c(i,j)A^ +1>j , (31) 

LgA™j = d(i, j)A™j-i + e{i,j)A™j + f(i,j)A™ j+1 . (32) 

To obtain a(i,j), b(i,j), c(i,j), d(i,j), e(i,j) and f(i,j), one 
has to use appropriate difference schemes for the operators 
L r and Lg . We use the Lax- Wendroff scheme for the ad- 
vective part, whereas the diffusive terms involving d 2 A/dr 2 
and d 2 A/d8 2 are handled by a simple Centred Space scheme 
(see, for example, Press et al. 1988). Since these diffusive 
terms are treated implicitly in one half-step and explicitly 
in the next half-step, what we get is essentially a Crank- 
Nicholson scheme for these diffusive terms. The terms in- 
volving dA/dr and dA/ 88 in the diffusive part are han- 
dled by an upwind scheme. The expressions of a(i,j), b(i,j), 
c(i,j), d(i,j), e(i,j) and f(i,j) are given below: 
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where U r (i,j) and Ug(i,j) respectively give the values of 
v r /r and vg / (r smO) at the grid point 

The equations are solved on a 64 x 64 grid. The bound- 
ary conditions are already discussed in the text. It was 
mentioned that the bottom boundary condition is to al- 
low the magnetic field to be carried freely through the bot- 
tom boundary, where the velocity is radially inward. We ad- 
vance A™j at the bottom by the simple upwind differencing 
scheme: 
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